function gr = gradp(FUN, x0 )
%  Purpose:    Computes the gradient vector or matrix (Jacobian) of a
%              vector-valued function that has been defined in a procedure.
%              Single-sided (forward difference) gradients are computed.
%
%  Format:     g = gradp( f, x0 ) ;
%
%  Input:      f    scalar, procedure pointer to a vector-valued function:
%
%                                          f:Kx1 -> Nx1
%
%              x0    Kx1 vector of points at which to compute gradient.
%
%  Output:     g     NxK matrix containing the gradients of f with respect
%                    to the variable x at x0.
%
%  Remarks:    gradp will return a row for every row that is returned by f.
%              For instance, if f returns a 1x1 result, then gradp will
%              return a 1xK row vector. This allows the same function to be
%              used where N is the number of rows in the result returned by f.
%              Thus, for instance, gradp can be used to compute the
%              Jacobian matrix of a set of equations.
%
%  Example:    proc myfunc(x);
%                 return_proc( x .* 2 .* exp( x .* x ./ 3 ) );
%              end_proc;
%
%              x0 = [ 2.5, 3.0, 3.5 ]' ;  is a 3x1 vector
%              y = gradp( myfunc, x0 );
%
%                           82.98901842    0.00000000    0.00000000
%                  y =       0.00000000  281.19752975    0.00000000
%                            0.00000000    0.00000000 1087.95414117
%
%              It is a 3x3 matrix because we are passing it 3 arguments and
%              myfunc returns 3 results when we do that.  The off-diagonals
%              are zeros because the cross-derivatives of 3 arguments are 0.
%
%  See Also:   hessp


% proc 1 = gradp(f,x0);
%     local f:proc;
%     local n,k,grdd,dh,ax0,xdh,arg,dax0,i,f0;

% /* check for complex input */
if ~isreal(x0)
    error('ERROR: Not implemented for complex matrices.');
end

f0 = feval( FUN, x0 ) ; 
n = rows(f0); % vector valued function FUN, Nx1.    %nota: rows and cols
k = rows(x0); % N.B.: column vector Kx1
grdd = zeros(n,k);

if cols( x0 ) > 1 %check x0 is a column vector
    error('ERROR: x0 must be a column vector.');
end

%/* Computation of stepsize (dh) for gradient */

ax0 = abs(x0);
if x0 ~= 0     % GAUSS: x0 /= 0;
    dax0 = x0./ax0;
else
    dax0 = ones(k,1);   %GAUSS = 1;
end

eps1 = 1e-8; % from Gauss
eps2 = 1e-2; % from Gauss
dh = eps1 * max( ax0,  eps2 * ones(k,1) ) .* dax0; %increments of x_0 Kx1 vector  
    %GAUSS 1e-8 * maxc( (ax0~(1e-2) * ones( rows(x0),1) )' ) .* dax0;
xdh = x0+dh; % Kx1
dh = xdh-x0;    % /* This increases precision slightly */
arg = diagrv( repmat(x0,1,k) , xdh ) ; % IO!
% GAUSS: arg = diagrv( reshape(x0,k,k)' , xdh ) ;

% note: RESHAPE takes columnwise elements from x0 and form a kxk matrix
% note: DIAGRV(A,v) puts the vector "v" in the main diagonal of the matrix "A"

i = 1;
while i <= k % GAUSS: do until i > k;
% for i = 1 : k, % oppure.
    grdd(:,i) = feval( FUN, arg(:,i) );
    i = i+1;
end

gr = ( grdd - repmat(f0,1,k) ) ./ repmat(dh',n,1) ; % forward derivative


